classdef NMRModel < handle
    %NMRMODEL Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        g_plant
        f_lowpass
        c_pid
        tList
        fList
        
        magneticNoise
        shotNoise
    end
    
    methods
        function obj = NMRModel(gamma, filter_order, filter_bw, varargin)
            %NMRMODEL Construct an instance of this class
            %   Detailed explanation goes here
                p=inputParser();
                p.addParameter('tList', 0:0.001:0.25, @isnumeric);
                p.addParameter('fList', logspace(-2, 3, 100), @isnumeric);
                p.parse(varargin{:});
                
                gain = ziBW2TC(gamma, 1)*360; % degree/Hz, 1s = 360 degree/Hz
                obj.g_plant = EquipmentControl.FeedBackOptimization.LTIUnit.LowPassFilter(1, gamma, 'gain', gain);
                obj.f_lowpass =EquipmentControl.FeedBackOptimization.LTIUnit.LowPassFilter(filter_order, filter_bw);
                
                obj.setTimeList(p.Results.tList);
                obj.setFrequencyList(p.Results.fList);
        end
        
        %% controller settings
        function obj = setTimeList(obj, tList)
            %METHOD1 Summary of this method goes here
            %   Detailed explanation goes here
            obj.tList = tList;
        end
        
        function obj = setFrequencyList(obj, fList)
            obj.fList = fList;
        end
        
        function obj = setPIDController(obj, pVal, iVal, dVal, bw)
            if nargin == 4
                Tf = 0;
            else
                Tf = ziBW2TC(bw, 1);
            end
            obj.c_pid = EquipmentControl.FeedBackOptimization.LTIUnit.PIDController(pVal, iVal, dVal, 'Tf', Tf);
        end
        
        function obj = setOptimizedController(obj, tList)
            if nargin == 2
                obj.setTimeList(tList);
            end
            
            tune = EquipmentControl.FeedBackOptimization.PIDOptimizer(obj.g_plant, obj.f_lowpass);
            tune.setControllerOptimal('type', 'PID', 'tList', obj.tList);
            [pVal, iVal, dVal] = tune.getPIDValue();
            obj.setPIDController(pVal, iVal, dVal);
        end
        
        %% noise settings
        function obj = setMagneticNoise(obj, a, bg, alpha)
            if nargin == 3
                alpha = 1; % default 1/f noise
            end
            
            obj.magneticNoise.alpha = alpha; % power index 
            obj.magneticNoise.a = a;         % strength of 1/f^alpha noise
            obj.magneticNoise.bg = bg;       % background white noise strength in [Hz/sqrt(Hz)]
            obj.magneticNoise.PSDfunc = @(f) a./(f.^alpha) + bg;
        end
        function psd = getMagneticNoise(obj, fList)
            if nargin == 1
                fList = obj.fList;
            end
            psd = obj.magneticNoise.PSDfunc(fList);            
        end
        
        function obj = setShotNoise(obj, bg)
            obj.shotNoise.bg = bg;  % phase shot noise in [deg/sqrt(Hz)]
            obj.shotNoise.PSDfunc = @(f) 0.0*f+bg;
        end
        function psd = getShotNoise(obj, fList)
            if nargin == 1
                fList = obj.fList;
            end
            psd = obj.shotNoise.PSDfunc(fList);
        end
        
        function [total_out_psd, phase_out_psd, magnetic_out_psd, bode_mag_Phase, bode_mag_Magnetic] = getFrequencyNoisePSD(obj, fList)
            if nargin == 1
                fList = obj.fList;
            end
            
            phase_in_psd = obj.getShotNoise(fList);
            magnetic_in_psd = obj.getMagneticNoise(fList);
            
            bode_mag_Phase = bode(obj.getPhaseNoise2FrequencyTransFunc, 2*pi*fList);
            bode_mag_Phase = bode_mag_Phase(:).';
            bode_mag_Magnetic = bode(obj.getReferenceTrackingTransFunc, 2*pi*fList);
            bode_mag_Magnetic = bode_mag_Magnetic(:).';

            
            phase_out_psd = phase_in_psd .* bode_mag_Phase;% .* additional_low_pass;
            magnetic_out_psd = magnetic_in_psd .* bode_mag_Magnetic;% .* additional_low_pass;
            total_out_psd = phase_out_psd+magnetic_out_psd;
        end
        
        function [mag, phase] = getSystemBode(obj, fList)
            if nargin == 1
                fList = obj.fList;
            end  
            [mag, phase] = bode(obj.getReferenceTrackingTransFunc, 2*pi*fList);
            mag = mag(:); phase = phase(:);
        end
        %% transfer functions
        function g = getPlantTransFunc(obj)
            g = obj.g_plant.getTFfunction();
        end
        function f = getFilterTransFunc(obj)
            f = obj.f_lowpass.getTFfunction();
        end
        function c = getControllerTransFunc(obj)
            c = obj.c_pid.getTFfunction();
        end
        
        function h = getReferenceTrackingTransFunc(obj)
            g = obj.getPlantTransFunc; f = obj.getFilterTransFunc; c = obj.getControllerTransFunc;
            h = feedback(g * f * c, 1);
        end
        function h = getPhaseNoise2FrequencyTransFunc(obj)
            g = obj.getPlantTransFunc; f = obj.getFilterTransFunc; c = obj.getControllerTransFunc;
            h = feedback(f * c, g);
        end
        function h = getNoise2PhaseTransFunc(obj)
            g = obj.getPlantTransFunc; f = obj.getFilterTransFunc; c = obj.getControllerTransFunc;
            h = feedback(f, g * c);
        end
        
        %% Response & Bode plot
        
    end
end

